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W. WALDMAN 
SUMMARY 

This technical memorandum describes the methods and computer programs used to specify, 
design, and implement time-domain and frequency-domain FIR digital filters for use in the 
analysis of F/A-18 flight test data. Two bandpass filters covering the 10-20 Hz and 32-52 Hz 
frequency bands have been developed, and they can be used for analysing the two dominant 
modes of structural vibration response occurring on the F/A -18 empennage. A highpass filter 
with a cutoff frequency of 8 Hz was also designed for filtering strain gauge data for use in 
producing fatigue load sequences for coupon testing. The effects of bandpass filtering on the 
transient response of short-term vibrations lasting about one second have also been 
investigated, leading to the conclusion that FIR digital filters have a negligible effect on the 
transient response characteristics of short bursts of vibration. 
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1.0 INTRODUCTION 


A joint effort between Canada and Australia to perform a full-scale structural fatigue test on the 
F/A-18 Hornet aircraft is currently underway. This project is called the International Follow On 
Structural Test Program (IFOSTP), and in 1989 a series of flight trials were performed at the 
Aerospace Engineering Test Establishment (AETE) in Canada to determine the Hornet vibration 
response under typical service conditions [1]. 

The airflow around the F/A-18 aircraft at moderate to high angles of attack is characterized by 
separated vortices produced by the highly-swept leading edge extensions (LEXes). When these 
vortices impinge on the horizontal and vertical tails, the buffet pressures produce high levels of 
structural vibration. In an attempt to reduce the levels of vibration response and the attendant 
fatigue damage, aerodynamic fences have been fitted to the leading edge extensions of all 
F/A-18 aircraft. This modification was designed by the manufacturer, McDonnell Douglas 
Corporation. 

The AETE flight trials measured typical sequences of dynamic loads on the vertical and 
horizontal tails, engine mounts, and other aft fuselage locations. The measurands consisted of 
accelerometer and strain gauge channels, and flight testing was carried out for ACM (Air 
Combat Manoeuvre) training, as well as for test points covering a range of AOA (angle of 
attack) and Q (dynamic pressure) conditions. The flights were conducted both with and without 
the LEX fence. 

The measurements that were collected are to be used to characterize the empennage vibration 
response environment experienced by the F/A-18. They will be processed to produce reference 
response spectra at a number of locations, and a suitable vibration loading is to be applied to the 
fatigue test article to match the in-flight vibration environment as closely as possible. Due to the 
broad-band nature of the vortex buffet pressure spectrum, there is significant structural response 
at the resonance frequencies of the vertical and horizontal tails. The 10-20 Hz and 32-52 Hz 
bands show a significant modal response, and the response in these two bands is commonly 
referred to as Mode I and Mode n, respectively. There is also a broad-band response in the 
80-100 Hz region under some AOA-Q conditions. The root-mean-square (RMS) response in a 
given band is a function of both AOA and Q, and the vibration levels of one mode relative to 
another also vary significantly. 

In order to be able to characterize the vibration response of the F/A-18, it is necessary to be able 
to determine the modal response in each of the important frequency bands. This calls for 
bandpass filtering of the data and, because it is stored in sampled serial digital form, the use of 
digital bandpass filters to perform the required signal processing is an efficient and convenient 
proposition. McDonnell Douglas performed flight tests during the development of the F/A-18 
aircraft, and their data analysis relied heavily on the use of bandpass filtering to determine the 
RMS response in the two main frequency regions of modal vibration. By filtering the AETE 
data into similar frequency bands, it will be possible to compare our results with existing data. 
Highpass filtering of strain gauge data is also required in order to leave only the dynamic 
response, which will later be combined with typical manoeuvre loads to produce fatigue load 
sequences for coupon testing. 

2.0 DESIGN OF BANDPASS AND HIGHPASS DIGITAL FILTERS 
2.1 Magnitude and phase response considerations 

Several parameters are used to describe a filter’s performance, and one of the most important of 
these is the magnitude response as a function of frequency. A filter's passband is the frequency 
region where the signal undergoes little or no modification, while the stopband is where the 
signal is attenuated to be less than some specified level. The transition band, as its name implies. 
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is the region between the passband and the stopband, and the signal in this region is attenuated 
to intermediate degree. When selecting a suitable filter, it is necessary to consider the flatness of 
response in the filter's passband, and limits on this are defined by how much variation or ripple 
we are prepared to tolerate. Another important parameter is the rate of attenuation in the 
transition band, which will determine how quickly a given amount of attenuation will be 
achieved. For many types of filters it is also necessary to consider the level of the ripple in the 
stopband. 

Evaluating the transfer function of a filter results in both a magnitude and a phase characteristic. 
Most filters have a non-linear phase shift characteristic, which causes each significant frequency 
component in the signal to undergo a different delay (non-linear group delay). This behaviour 
introduces distortion of the output waveform in addition to the simple removal of unwanted 
frequency components. As a result, the peaks in the output waveform may be different in the 
presence of delay distortion than without it, and this can affect the results of peak-valley 
processing as commonly carried out in fatigue damage studies. 

2.2 FIR digital filters 

In terms of magnitude and phase response, it is possible to design digital filters that possess the 
same characteristics as classical analogue filter designs. For data that have already been 
measured and are stored on computer disk, it is possible to completely eliminate the non-linear 
phase behaviour by processing the data in both the forward and reverse directions. This 
procedure exactly cancels the phase shifts, and also results in a squaring of the magnitude 
response function. 

There exists a class of digital filters that can be designed to exhibit a linear phase behaviour, 
together with very high rates of attenuation in the transition band. The latter property is 
particularly desirable to enable the separation of F/A-18 vibration response data into its 
individual modal responses with as little interaction between bands as possible. These filters are 
called finite impulse-response (FIR) linear phase filters, and a Fortran program which enables 
optimum equi-ripple FIR lowpass, highpass, bandpass and bandreject filters to be designed is 
available and is based on the McClellan-Parks-Rabiner algorithm [2]. An updated version of this 
program appears in [3], and is the version that was eventually adopted for the FIR filter design 
work. 

In order to enable the computed filter coefficients to be easily verified, the program described in 
[3] was modified to calculate and display the magnitude and phase of the Fourier transform of 
the filter impulse response. The group delay is also computed and output, making it easy to 
check the filter's linear phase performance in the passband. The program implemented here is 
called EQFIR. 

Once the filter design has been completed and a set of suitable filter coefficients determined, the 
filtering action can be obtained by time-domain convolution of the input signal with the filter's 
symmetrical impulse response. An efficient algorithm for performing time domain filtering has 
been presented by Rabiner [4], together with Fortran source code for implementation as a 
callable function. This function was modified to make it possible to apply it to simultaneously 
filtering multiplexed channels of data, as found in AETE data files. The new version is 
presented in Appendix A and the function is called FILT1. An example program illustrating the 
use of this function is presented later. 

2.3 Design specification for FIR digital filters 

In order to design an FIR filter using program EQFIR, it is necessary to prepare a specification 
that defines the filter performance characteristics. Figure 1 shows a general bandpass filter that 
can be used as a basis for developing a suitable design specification. M re is the desired 
magnitude response in the passband, and is usually set to unity. The passband is specified by the 
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Figure 1: Parameters for FIR filter specification 


two edge frequencies of the interval [Q P1 , QpJ, and the two stopbands are specified by the edge 
frequencies for the intervals [0, £2 S iI and [Q^, 0.50]. Note that the frequency 12 = co/co s has been 
normalized by the sampling frequency, co s = 27tf s , so 12 = 0.50 corresponds to the Nyquist 
frequency. The desired ripple in the passband is denoted by 8 j, and is measured as a peak value 
relative to M ra . The ripple in the stopband is denoted by 83 , and is measured as a peak value 
relative to the zero level. 

If multiple transition bands are present in the filter, it is important to note that they must all be of 
the same width. If they are not, then the filter design program may not be able to achieve a 
smooth rolloff through the transition bands, resulting in unsatisfactory filter characteristics. 

2.4 Design procedure for FIR digital filters 

The FIR filter design program EQFIR expects four lines of data to be provided, and the required 
format is presented below: 

NFILT JTYPE NBANDS LGRID SFREQ 
EDGE(1) ... EDGE(2*NBANDS) 

DESR(l) ... DESR(NBANDS) 

WGHT(l) ... WGHT(NBANDS) 

Each of the parameters required by the specification of the filter design problem are defined as 
follows: 

NFILT The filter length in samples. NFILT must satisfy 3 £ NFILT S NFMAX. In 
the present implementation NFMAX has been set to 1024. 
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JTYPE The type of filter. 


JTYPE = 1 for multiple passband/stopband. 

JTYPE = 2 for differentiator. 

JTYPE = 3 for Hilbert transformer. 

The last two filter types are outside the scope of this memorandum and will 
not be discussed any further. 

NBANDS The number of frequency bands, up to a maximum of 10 bands. 

LGRID The grid density in each band used for computing the extremal frequencies, 
assumed to be 16 if a value of zero is specified. 

SFREQ The sampling frequency for the filter design. Used to scale the frequency 
axis for listings of filter frequency response and group delay characteristics. 

EDGE An array of size 20, containing the frequency bands, specified by upper and 
lower cutoff frequencies, up to a maximum of 10 bands. 

DESR An array of size 10, containing the desired frequency response in each 
band. In the passband the desired value is usually set to unity, and in the 
stopband it is usually set to zero. 

WGHT An array of size 10, containing the positive weight function for each band. 

Each value is the ratio of the passband ripple and ripple in the given band. 
For stopbands it is usually calculated as 8 ,/ 8 2 , while for passbands it is 
usually set to unity. 

Note that the frequencies to be input in the EDGE array are normalized by the sampling 
frequency, and each frequency band pair specifies a closed subinterval of the frequency axis 
[0, 1/2]. The arrays DESR and WGHT then specify the ideal desired response and the weight 
function in each band. 

When designing a digital filter using EQFIR, it is necessary to specify the length of the filter in 
samples. A program for computing an estimate for the value of N is presented in Appendix B. It 
is based on a formula developed for estimating the value of N for lowpass digital FIR filters [5], 
and can also be applied to estimating N for bandpass filters. Given the filter parameters Q s , flp, 
5, and 83 , the minimum filter impulse response duration required to meet the specifications can 
be estimated from the relation 

N = D^f^SjVAQ - (0.51244 log , 0 K + 11.01217)Afl + l 

where 

K = 8,/82 

Doo( 8 „ 62 ) = [ 0.005309 (log 10 8 ,)* + 0.07114 (log , 0 8 ,) - 0.4761 ] log I0 8 2 - 
( 0.00266 (log I0 83)2 + 0.5941 (log 10 82 ) + 0.4278 ] 

A£2= iQp-OsI 

The above equation generally provides a good first estimate of N, particularly for low-order FIR 
filters, as the relative error is 1.3% for 8 | £0.1 and S 2 £0.1 [5]. If the iterated filter design 
obtained using EQFIR does not meet the required filter specification, then the value of N can be 
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increased and the filter coefficients re-computed. Using this technique of successive 
approximation, a satisfactory solution is usually obtained after only a few additional trials. 

The theoretical group delay G of an FIR filter can be quickly deter min ed by using the following 
formula 

(N- 1) 

G = - (seconds) 

2 f s - 


When designing analogue filters, it is common to specify the rolloff rate in the transition band in 
terms of dB/octave. For n-pole Butterworth filters, the rolloff rate is simply 6 ndB/octave 
(20n dB/decade), which yields a linear plot on a logarithmic frequency scale. When designing 
digital filters using EQFIR, it is the width of the transition band that has to be specified. 

If the required rate of attenuation in dB/octave is known, and the magnitude response in the 
passband is assumed to be unity, then the width of the lower and upper transition bands, Aft LTB 
and AfluTB, can be estimated from the following two formulae 

AQltb = QpiU + R-ltb/( 20 log 10 Sj - Rltb)) 

A^utb = ~ Gpj(20 log 10 ^/Rutb 

where 

Rltb = desired average rolloff rate in the lower transition band in dB/octave 

Ruts = desired average rolloff rate in the upper transition band in dB/octave 

For a given rolloff rate, the width of the lower transition band will always be larger than the 
width of the upper transition band. If a minimum average rolloff rate is required, then this will 
enforce a constraint on the maximum width of the lower transition band. The edge frequencies 
of the upper and lower transition bands can be computed using the relations 

Qsi = — A£2 L tb 

Os2 = ^P2 + 


It is also common practice in filter design to specify the passband and stopband ripple in dB 
relative to the passband amplitude. Hence it is necessary to be able to convert such 
specifications into a format compatible with the input parameters of EQFIR. If 5 dB1 and 8 dB2 are 
the passband and stopband ripple specified in dB relative to the passband amplitude, then the 
corresponding values of 8 ! and 82 can be computed from the following relations 


SdBi/20 

8 ,= 10 


1 


$dB2/20 

8j= 10 


2.5 Design of a 10-20 Hz bandpass FIR filter 

In order to isolate the Mode I response in the 10-20 Hz band, a digital FIR bandpass filter was 
designed using the techniques described above. The sampling rate was f s = 606.06 Hz, and the 
bandpass edge frequencies were taken to be 10 Hz and 20 Hz. The passband ripple was specified 
to be no greater than 0.05 dB, while the stopband ripple was to be -60 dB or better. In order to 
provide good rejection of any low-frequency response, the desired average rolloff rate in the 
lower transition band was specified to be 140 dB/octave. Hence we have that 
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f $ = 606.06 Hz 


8 <ibi ~ 

0.05 dB 

<=> 

5, = 0.00577 

SdB2 = 

-60 dB 

<=> 

83 = 0.001 

Rltb = 

140dB/octave 



II 

.5! 

10 Hz 

<=> 

a M = 0.01650 

fp2 = 

20 Hz 

<=> 

Q K = 0.03300 

A^ltb - 

. 3 Hz 

<=> 

AQ ltb = 0.00495 

fsi = 

7 Hz 

<=> 

Qs,= 0.01155 

*■*» 

3 

11 

23 Hz 

<=> 

Qs 2 = 0.03795 


Now, in the formula for estimating N, AI2 = A£2 LTB = 0.00495, and using the above values of S, 
and 62 we obtain 

N = 549 

K = 5.77 

The data file for EQFIR for the 10-20 Hz bandpass filter then becomes 


549 1 

3 

16 

606.06 


0.0 0.01155 

0.01650 0.03300 

0.03795 0 

.5 

0.0 1.0 

0.0 




5.77 1.0 

5.77 




results output by program EQFIR for N = 549 are presented below. 




BAND 1 

BAND 2 

BAND 3 

LOWER BAND EDGE 

HZ 

0.000 

10.000 

23.000 

UPPER BAND EDGE 

HZ 

7.000 

20.000 

303.030 

LOWER BAND EDGE 


0.0000000 

0.0165000 

0.0379500 

UPPER BAND EDGE 


0.0115500 

0.0330000 

0.5000000 

DESIRED VALUE 


0.0000000 

1.0000000 

0.0000000 

WEIGHTING 


5.7700000 

1.0000000 

5.7700000 

DEVIATION 


0.0014360 

0.0082859 

0.0014360 

DEVIATION IN DB 


-56.8567657 

0.0716737 

-56.8567657 


Looking at the bottom two rows giving the deviation, it can be seen that this filter does not quite 
meet the design specifications, as both the passband and stopband ripple are slightly greater than 
desired. After a number of additional trials a final value of N = 571 was obtained, and the results 
of the design process are given below. 


LOWER BAND EDGE HZ 
UPPER BAND EDGE HZ 
LOWER BAND EDGE 
UPPER BAND EDGE 
DESIRED VALUE 
WEIGHTING 
DEVIATION 
DEVIATION IN DB 


BAND 1 
0.000 
7.000 
0.0000000 
0.0115500 
0.0000000 
5.7700000 
0.0008924 
-60.9890862 


BAND 2 
10.000 
20.000 
0.0165000 
0.0330000 
1.0000000 
1.0000000 
0.0051490 
0.0446090 


BAND 3 
23.000 
303.030 
0.0379500 
0.5000000 
0.0000000 
5.7700000 
0.0008924 
-60.9890862 


The upper plot in Figure 2 shows the symmetric impulse response of the filter. The middle plot 
shows the magnitude response of the filter over a frequency range of 0-100 Hz, and it is seen 
that the passband ripple is within 0.05 dB, and the stopband ripple is less than -60 dB. The 
lower plot shows the group delay characteristics of the bandpass filter computed from the filter’s 
phase response, and it is seen that a uniform group delay of 0.4703 seconds is obtained in the 
passband. 
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Figure 2: Impulse response, magnitude response and group delay of 10-20 Hz FIR 
bandpass filter. 
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The time delay indicated by the lower plot corresponds closely to the group delay 
G = 0.4711 seconds computed from the theoretical formula presented earlier. The group delay 
plot confirms that the filter exhibits linear-phase characteristics in the passband. This group 
delay is also evident in the stopband, except that a 180* phase change is present at each of the 
zeros in the transfer function. 


2.6 Design of a 32-52 Hz bandpass FIR filter 

In order to isolate the Mode II response in the 32-52 Hz band, a digital FIR bandpass filter was 
designed to work with the data sampled at f s = 606.06 Hz, and the bandpass edge frequencies 
were taken to be 32 Hz and 52 Hz. The passband ripple was specified to be no greater than 
0.05 dB, while the stopband ripple was to be no greater than -60 dB. In order to provide good 
rejection of any low-frequency response, particularly any Mode I response that might be present, 
the desired average rolloff rate in the lower transition band was specified to be 144dB/octave. 
Hence we have that 


f s = 606.06 Hz 


8 dB i= 0.05 dB 


8 > = 0.00577 

5dB2 = ~60 dB 

<=> 

82 = 0.001 

Rltb - 14^ dB/octave 



f P1 = 32 Hz 

«=> 

Q P1 = 0.05280 

f M = 52 Hz 

<=> 

ftp 2 = 0.08580 

Af LTB = 9.41Hz 

<=> 

An LTB = 0.01553 

f sl = 22.59 Hz 

<=> 

Qs, = 0.03727 

f s2 = 61.41Hz 

<=> 

Qs 2 = 0.10133 


Now, AQ = AQ ltb = 0.01553, and using 8 j and 82 in the formula for estimating N we obtain 
N = 176 
K = 5.77 

A data file suitable for use with program EQFIR was created, and the resulting filter for N = 176 
did not quite meet the design specifications for stopband and passband ripple. After additional 
trials with increased values of N, a suitable design was obtained with N = 195. The final version 
of the data file is given below. 


195 1 

3 

16 

606.06 


0.0 0.03727 

0.05280 0.08580 

0.10133 0.5 


0.0 1.0 

0.0 




5.77 1.0 

5.77 




results obtained by program EQFIR for N = 195 are presented below. 




BAND 1 

BAND 2 

BAND 3 

LOWER BAND EDGE 

HZ 

0.000 

32.000 

61.412 

UPPER BAND EDGE 

HZ 

22.588 

52.000 

303.030 

LOWER BAND EDGE 


0.0000000 

0.0528000 

0.1013300 

UPPER F AND EDGE 


0.0372700 

0.0858000 

0.5000000 

DESIRED VALUE 


0.0000000 

1.0000000 

0.0000000 

WEIGHTING 


5.7700000 

1.0000000 

5.7700000 

DEVIATION 


0.0009579 

0.0055273 

0.0009579 

DEVIATION IN DB 


-60.3732834 

0.0478770 

-60.3732834 


The upper plot in Figure 3 shows the symmetric impulse response of the filter. The middle plot 
shows the magnitude response of the filter over a frequency range of 0-100 Hz. and it is seer, 
that the passband ripple is within 0.05 dB, and the stopband ripple is less than -60 dB. 
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Figure 3: Impulse response, magnitude response and group delay of 32-52 Hz FIR 
bandpass filter. 
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The lower plot shows the group delay characteristics computed from the bandpass filter's phase 
response, and it is seen that a uniform group delay of 0.1601 seconds is obtained in the passband 
which, to all intents and purposes, is identical to the theoretical group delay. 


2.7 Design of an 8 Hz cutoff highpass FIR filter 

In order to filter the measured vibration response to remove any very low-frequency or static 
response, a digital FIR highpass filter with 8 Hz cutoff was designed to work with the data 
sampled at f s = 606.06 Hz. The cutoff frequency was chosen so that there would be no 
significant filtering of any Mode I response present in the measured data. The passband ripple 
was specified to be no greater than 0.05 dB, while the stopband ripple was to be no greater than 
-60 dB. In order to provide good rejection of any low-frequency response, the upper edge of the 
stopband was set to 5 Hz. 


f s = 606.06 Hz 



5^,= 0.05 dB 

<=> 

6 , = 0.00577 

5 d B2 — -60 dB 

t=> 

82 = 0.001 

f sl = 0 Hz 

<=> 

fts, = 0.00000 

f S2 = 5 Hz 

<=> 

fts 2 = 0.00825 

f P i = 8 Hz 

<=> 

ftp, = 0.01320 

fp 2 = 303.03 Hz 

<=> 

ftp 2 = 0.50000 

Af-re = 3 Hz 


Aft TB = 0.00495 


Now, Aft = Aftre = 0.00495, and using 8 , and 82 in the formula for estimating N we obtain 
N = 549 
K - 5.77 


A data file suitable for use with program EQFIR was created, and the resulting filter for N = 549 
did not quite meet the design specifications for stopband and passband ripple. After additional 
trials with increased values of N, a suitable design was obtained with N = 559, and the final 
version of the data file is given below. 


559 

1 

2 

0.0 

0.00825 

0 . 

0.0 

1.0 

0 . 

5.77 

1.0 

5. 


16 606.06 
01320 0.5 

0 

77 


The results output by program EQFIR for N = 559 are presented below. 


LOWER BAND EDGE HZ 
UPPER BAND EDGE HZ 
LOWER BAND EDGE 
UPPER BAND EDGE 
DESIRED VALUE 
WEIGHTING 
DEVIATION 
DEVIATION IN DB 


BAND 1 
0.000 
5.000 
0.0000000 
0.0082500 
0.0000000 
5.7700000 
0.0009408 
-60.5296059 


BAND 2 
8.000 
303.030 
0.0132000 
0.5000000 
1.0000000 
1.0000000 
0.0054287 
0.0470253 


The upper plot in Figure 4 shows the symmetric impulse response of the filter. The middle plot 
shows the magnitude response of the filter over a frequency range of 0-100 Hz, and it is seen 
that the passband ripple is within 0.05 dB, and the stopband ripple is less than -60 dB. The 
lower plot shows the group delay characteristics computed from the highpass filter’s phase 
response, and it is seen that a uniform group delay of 0.4604 seconds is obtained in the 
passband. This is in close agreement with the theoretical group delay G = 0.4579 seconds. 
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Figure 4: Impulse response, magnitude response and group delay of 8 Hz cutoff FIR 
highpass filter. 
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Figure 5: Bandpass filtered responses of multi-frequency signal. 


3.0 FILTERING PERFORMANCE OF FIR FILTERS 

Convolution of a filter's impulse response with the input signal in order to perform the filtering 
operation is a process that is very computationally intensive. This belies the relative simplicity 
of the time-domain FIR filtering algorithm. In order to reduce the computation time, a frequency 
domain technique using FFTs has been developed [6]. For larger filter lengths it offers 
significant speed gains, which is beneficial when filtering large amounts of flight test data. 

A Fortran subroutine implementing the "overlap-add" technique of fast convolution in the 
frequency domain can be found in [3]. The name of the subroutine is RFILT, and a listing is 
provided in Appendix C. 

A Fortran program TESTFELT was written in order to verify the operation of the filtering 
algorithms when using the filter coefficients produced by program EQFIR, and a listing is 
provided in Appendix D. TESTFILT generates a multi-frequency signal consisting of four equal 
amplitude sine waves with frequencies at 6 Hz, 15 Hz, 45 Hz and 90 Hz. These frequencies were 
chosen to be representative of those that occur in F/A-18 empennage vibrations. Both time- 
domain and frequency-domain filtering techniques were used to permit a comparison of filtering 
quality, as well as providing an independent check of the results of each method. 

The upper plot in Figure 5 shows the waveform produced using the four equal amplitude sine 
waves. The middle and lower plots show the results after applying the 10-20 Hz and 32-52 Hz 
bandpass filters. Note that the results of filtering in the time and frequency domains were 
effectively identical to about five significant figures, so there was no discernible difference 
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Figure 6: Comparison of simulated 15 Hz response transient with the 10-20 Hz 
bandpass filtered response. 

between the results when the plots were overlaid. The bandpass filtered responses also show 
clearly the time delay introduced by the filtering process. 

For the 10-20 Hz bandpass filter, it was found that the frequency-domain filtering technique 
was 7.7 times faster than its time-domain counterpart. The frequency-domain 32-52 Hz 
bandpass filter was found to be 2.4 times faster than the time-domain implementation. For a 
constant size of FFT block length, the speed of the frequency-domain filter was independent of 
the number of filter coefficients, while the speed of the time domain filter was approximately 
inversely proportional to the number of coefficients. 

The frequency-domain filtering technique relies on the swiftness of the forward and inverse 
FFTs to obtain its speed improvement. A number of other implementations of the FFT algorithm 
were incorporated into subroutine RFILT, but the original code using subroutines FAST and 
FSST (given in [3]) proved to be the most efficient. As a result of their computational efficiency 
and well-documented performance characteristics, these two subroutines were also adopted for 
all other work involving FFT computations. 

4.0 TRANSIENT PERFORMANCE OF FIR RLTERS 

Filtering operations are normally accompanied by changes in the waveform of the signal, 
particularly if the signal has many frequency components and some of these are removed or 
significantly attenuated. Before applying the FIR filters designed above to F/A-18 vibration 
responses, the effects of the filtering action on a simulated pulse-type transient waveform that 
represents a burst of vibration at a particular modal frequency was studied. The Fortran program 
CHECKQ in Appendix E can generate a 15 Hz or 45 Hz sinusoidal waveform that is 
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Frequency in Hz 


Figure 7: Comparison of power spectra for simulated 15 Hz response transient before 
and after 10-20 Hz bandpass filtering. 


approximately one second in length. To make this pulse better simulate the build-up and decay 
of a typical structural vibration response transient, the waveform can be windowed using an 
exponential or a sine-squared weighting function. The damping ratio of the decay in the 
exponential window was chosen to be 3% of critical, which is typical of aircraft structural 
vibration modes. The program then filters the transient, and the output consists of listings of the 
time histories and power spectra of the unfiltered and filtered responses. The time history listing 
includes the filtered sequence with the theoretical time delay removed, making it easy to 
conduct a point by point comparison of the filtered waveform with the unfiltered data. 

The upper plot in Figure 6 shows the unfiltered response of the simulated 15 Hz exponentially- 
weighted response transient, and the middle plot shows the resulting 10-20 Hz FIR bandpass 
filtered response. The lower plot shows the time-delay corrected version of the filtered response 
superimposed over the unfiltered transient. It is evident that the bandpass filtered response is 
very similar to the raw data, except at the leading and trailing edges of the pulse where there is a 
little extra oscillation that is not present in the unfiltered data. This is caused by the removal of 
some of the frequency components associated with the discontinuity at the leading and trailing 
edges of the simulated pulse. 

Figure 7 shows the power spectra obtained from the time histories of the unfiltered and bandpass 
filtered 15 Hz pulse. In the 10-20 Hz passband of the bandpass filter, there is very good 
agreement between the two spectra, with all the side lobes being present with the correct 
amplitude. As expected, the energy in the stopbands has been significantly reduced by the action 
of the bandpass filter. 

The results shown in Figures 6 and 7 show that the processing performed by the FIR filters will 
provide only a minor modification to any transient waveforms falling within the passband of the 
filter. Even better results were obtained for the sine-squared weighting of the sinusoidal pulse, 
but these are not presented here as it was found that the sine-squared weighting had created a 
less demanding transient because it goes smoothly to zero at either end. 
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5.0 CONCLUSION 


A complete design procedure for computing FIR digital filter coefficients has been presented. A 
set of filters suitable for processing F/A-18 flight test data has been determined, and includes 
10-20 Hz and 32-52 Hz bandpass filters, as well as an 8 Hz highpass filter. A program for 
performing time-domain and frequency-domain filtering was written and has been used to verify 
the filter coefficients and the operation of the filtering algorithms. A program for studying the 
effects of bandpass filtering on a simulated pulse-like sinusoidal transient waveform whose 
fundamental frequency lies in the filter's passband was also written. The results indicate that the 
FIR digital filters used here have a minimal effect on the waveforms of transients whose 
frequency content lies within the passband of the filter. 
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APPENDIX A 


A Fortran Subroutine For Digital Filtering In The Time Domain 


The following Fortran subroutine can be used to perform efficient digital filtering of signals in 
the time domain. The coefficients of the filter impulse response can be computed using program 
EQFIR. This subroutine is based on Rabiner’s algorithm [4], but has been modified to make it 
easy to apply to multiple channels of data contained in the one data file. 


A subroutine that can read in filter coefficients from an output file written by the filter design 
program EQFIR is also provided, and can be found at the end of the listing. 

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 


c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


REAL FUNCTION FILT1(X,XSAV,H,N,IPT,INIT) 

X A SINGLE, REAL DATA SAMPLE 

XSAV A WORK ARRAY, DIMENSION OF 2*N 

H FILTER COEFFICIENT ARRAY, DIMENSION OF N 

N NUMBER OF FILTER COEFFICIENTS 

IPT A LOCAL VARIABLE THAT MUST BE SAVED BETWEEN 

SUCCESSIVE CALLS TO THIS FILTER ROUTINE. 

A SEPARATE ARGUMENT WAS ADDED TO ALLOW THIS 
ROUTINE TO BE USED WHEN FILTERING MORE THAN 
ONE SEQUENCE AT A TIME. 

INIT=0 TO INITIALIZE FILTER BEFORE RUNNING 
-1 WHEN RUNNING 

FILT1 ON RETURN, FILT1 IS THE FILTERED OUTPUT SAMPLE 

REAL X,H(1),XSAV(1),Y 
INTEGER N,INIT,N2,IPT,I 

IF (INIT.EQ.O) THEN 
DO I * 1,2 *N 
XSAV(I) =0.0 
ENDDO 

IPT = N + 1 
FILT1 =0.0 
RETURN 
ENDIF 


XSAV(IPT) = X 
XSAV(IPT-N) = X 
Y-0.0 

DO I = 1,N 

Y = Y + H(I)*XSAV(IPT-I+1) 
ENDDO 

IPT = IPT + 1 

IF (IPT.GT.2*N) IPT = N + 1 
FILT1 = Y 


RETURN 

END 


CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 

SUBROUTINE READH(NFILT,H,FILENAME,UNITNO) 

C THIS PROCEDURE READS IN THE COEFFICIENTS OF THE FILTER IMPULSE 
C RESPONSE INTO THE ARRAY H. 


16 






DIMENSION H (*) 

INTEGER NFILT,UNITNO,IPOS 

CHARACTER FILENAME*(*) 


OPEN(UNIT-UNITNO,FILE-FILENAME,STATUS- 1 OLD’,READONLY) 

CALL FINDSTR(UNITNO,'FILTER LENGTH - ’,IPOS) 

READ(UNITNO,'(35X,15)') NFILT ! READ FILTER LENGTH 

DO I - 1,3 ! SKIP 3 LINES 

READ(UNITNO,*) 

ENDDO 

! READ FILTER COEFFICIENTS 
DO I - 1, (NFILT+l)/2 

READ(UNITNO,'(15X,14,4X,E15.8,5X,14)') ITEMP1,RTEMP1,ITEMP2 
H(ITEMPl) - RTEMP1 
H(ITEMP2) - RTEMP1 
ENDDO 

CLOSE(UNIT—UNITNO) 

RETURN 

END 


CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 

SUBROUTINE FINDSTR(UNITNO,SUBSTR,IPOS) 

INTEGER UNITNO,IPOS 
CHARACTER SUBSTR*(*) 

CHARACTER LINE*132 
INTEGER INDEX 

IPOS - 0 

DO WHILE (IPOS.EQ.O) 

READ(UNITNO,'(A132)') LINE 
IPOS - INDEX(LINE,SUBSTR) 

ENDDO 

IF (IPOS.GT.O) BACKSPACE(UNITNO) 

RETURN 

END 
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APPENDIX B 


A Fortran Program For Estimating N For FIR Digital Filters 


The following Fortran program can be used for estimating the value of N for FIR digital filters. 
The required input data consists of the sampling rate f s , the width of the transition band Af-^, 
and the passband and stopband ripple 5] and 65 . Note that the frequency information can be 
supplied in any consistent units, as it is automatically non-dimensionalised with respect to f s 
within the program. 

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

PROGRAM ESTIMN 

C This program is used to estimate the order of the FIR filter 
C required to meet the filter specifications. 

IMPLICIT NONE 

INTEGER N 

REAL FSAMPLE,DELTAF,DELTAl,DELTA2 

REAL A1,A2,A3,A4,A5,A6,B1,B2,DF, F,D,LOGD1,LOGD2 

REAL ALOG10 

WRITE(*,'(/,IX,A,$)’) 'INPUT Fs, Delta_F, Delta_l andDelta_2: ' 

READ(*,*) FSAMPLE,DELTAF,DELTAl, DELTA2 

A1 - 5.309E-3 

A2 - 7.114E-2 

A3 - -4.761E-1 
A4 - -2.660E-3 
A5 - -5.941E-1 
A6 - -4.278E-1 

B1 - 11.01217 
B2 - 0.51244 

DF - DELTAF/FSAMPLE 

LOGD1 - ALOGIO(DELTAl) 

LOGD2 - ALOGIO(DELTA2) 

F - B1 + B2MLOGD1 - LOGD2) 

D - (A1*LOGD1*LOGD1 + A2*LOGDl + A3)*LOGD2 

& + A4*LOGD1*LOGD1 + A5*LOGDl + A6 

N - D/DF - F*DF + 1 

WRITE(*,’(/,IX,A,15,/)') ’Estimated N = ',N 

STOP 

END 
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APPENDIX C 


A Fortran Subroutine For Digital Filtering In The Frequency Domain 


The Fortran subroutine RFILT [3] can be used to perform digital filtering of signals in the 
frequency domain. A listing of this subroutine is presented below. 

SUBROUTINE: RFILT 

FILTER ONE FRAME (I.E., NP POINTS) OF DATA 
PROGRAM ASSUMES: 

1. R(NP+1) TO R(2*NP) HAS NOT BEEN USED SINCE LAST CALL TO 
'RFILT' AND WAS ZERO BEFORE FIRST CALL. 

2. INPUT AND OUTPUT DATA (TIME SERIES) ARE IN R(l) TO R(NP) 

3. F(1) TO F(2*NP+2) IS THE FILTER FREQUENCY RESPONSE. 

4. DATA ARE STORED IN F AS: DC,0.,F1REAL,F1IMAG, 

F2REAL, F2IIIAG, . . ., FSAMPLING/2 ., 0 . 

5. S(l) TO S(NP) IS A SCRATCH ARRAY; IT MAY BE USED IN 'MAIN'. 

6. IMPULSE RESP. OF FILTER F IS ZERO FOR i:P+l THRU 2*NP 
(I.E. SECOND HALF OF H(T)-0, WHERE F(W)-FFT(H) ) 


SUBROUTINE RFILT(R, F, S, NP) 

DIMENSION S(l), R(1), F(l) 

NPT2 - NP*2 

STORE PREVIOUS TAIL IN SCRATCH ARRAY S(.); 

ZERO SECOND HALF OF R—FIRST HALF OF R IS NEW DATA TO BE FILTERED 

DO 10 K—1,NP 
NPPK - NP + K 
S(K) - R(NPPK) 

R(NPPK) - 0. 

10 CONTINUE 

TAKE FFT OF DATA 

CALL FAST(R, NPT2) 

HANDLE VALUES AS COMPLEX VALUES 

KMAX - NPT2 + 1 
DO 20 K-1,KMAX,2 

X - F(K)*R(K) - F(K+l)*R(K+l) 

Y - F(K)*R(K+l) + F(K+l)*R(K) 

R(K) - X 
R(K+l) - Y 

20 CONTINUE 

INVERSE TRANSFORM PRODUCT 
CALL FSST(R, NPT2) 

ADD IN TAIL FROM PREVIOUS FRAME WHICH WAS STORED IN S(.) 

DO 30 K«1,NP 

R(K) - R(K) + S(K) 

30 CONTINUE 
RETURN 
END 
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APPENDIX D 


A Fortran Program For Testing FIR Filter Designs 


The Fortran program TESTFILT can be used to test the efficacy of the FIR filtering algorithms 
using filter coefficients that have been generated by the filter design program EQFIR. The 
listing shown below is an example of how to implement digital filtering in the time-domain and 
frequency-domain using the Fortran function and subroutine presented in Appendices A and C. 

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

PROGRAM TESTFILT 

C This program filters a time series point by point in the time 
C domain and in the frequency domain. 

C 

C The execution times for both types of filtering are displayed 
C for the purpose of easy comparison. 

IMPLICIT NONE 

INTEGER NH,NHT2,NFAST,NDATA 
REAL SRATE,PI,TWOPI,DT 


PARAMETER ( NH = 1024 ) 
PARAMETER ( NHT2 = 2*NH ) 
PARAMETER ( NFAST = 2*NHT2+2 ) 
PARAMETER ( NDATA = 3*NHT2 ) 
PARAMETER ( SRATE « 606.06 ) 
PARAMETER ( DT = 1.0/SRATE ) 
PARAMETER ( PI - 3.141592654 ) 
PARAMETER ( TWOPI = 2.0*PI ) 


REAL* 8 FSB1(10),FPB<10), FSB2(10),T 

REAL XSB1,XSB2,XPB,WT,CPUTD,CPUFD 

REAL H(NH),XSAV(NHT2),F(NFAST),R(NFAST),WORKl(NHT2) 

REAL DTOT(NDATA),DPB(NDATA),DFIL(NDATA),DFILT(NDATA) 

INTEGER NFILT,IFILTER,I,J,NSB1,NSB2,NPB 
INTEGER IB r NBLOCKS,IPART,IPOS r NPT 
CHARACTER FILENAME*20 

! EXTERNAL FUNCTIONS 

REAL SIN,SECNDS,FILT1 

REAL*8 DMOD 

INTEGER JMOD 

IFILTER - 0 

DO WHILE (IFILTER.EQ.0) 

WRITE {*,*) 

WRITE(*,*) ’PROGRAM FOR TESTING FIR FILTER DESIGNS' 

WRITE(*,*) *---' 

WRITE(*,*) ' 1 * 10 HZ TO 20 HZ BANDPASS FILTER' 

WRITE(*,*) ' 2 » 32 HZ TO 52 HZ BANDPASS FILTER' 

WRITE(*,*) ' 3 - 5 HZ LOWPASS FILTER’ 

WRITE(*,*) ' 4 - 8 HZ HIGHPASS FILTER' 

WRITE(*,*) 

WRITE (*, ' (1X,A,$) ') 'INPUT YOUR SELECTION: ' 

READ(*,*) IFILTER 

IF (IFILTER.LT.l .OR. IFILTER.GT.4) IFILTER = 0 
ENDDO 

IF (IFILTER.EQ.1) THEN 
FILENAME - 'BP15.0UT' 

NSB1 - 1 
NPB - 1 
NSB2 » 2 
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FSBl(l) - 6.0 

FPB(l) - 15.0 

FSB2(1) - 45.0 

FSB2(2) - 90.0 

ELSE IF (IFILTER.EQ.2) THEN 


FILENAME - 

'BP45 

NSB1 - 2 


NPB - 1 


NSB2 - 1 


FSBl(l) - 

6.0 

FSB1(21 - 

15.0 

FPB(l) - 

45.0 

FSB2(1) - 

90.0 


ELSE IF (IFILTER.EQ.3) THEN 
FILENAME - 'LP05.OUT' 

NPB - 1 
NSB1 - 1 
NSB2 - 0 
FPB(l) - 2.0 

FSBl(l) - 200.0 
ELSE IF (IFILTER.EQ.4) THEN 
FILENAME - ’HPOS.OUT' 

NPB - 1 
NSB1 - 1 
NSB2 - 0 
FPB(l) - 15.0 

FSBl(l) - 5.0 

END IF 

! CREATE DATA SEQUENCE 

DO I - 1,NDATA 
T - (1-1)*DT 
XSB1 - 0.0 
XPB -0.0 
XSB2 - 0.0 
DO J - 1,NSB1 

WT - TWOPI*DMOD(FSBl(J)*T,1.0D0) 

XSB1 - XSB1 + SIN(WT) 

ENDDO 

DO J - 1,NPB 

WT - TWOPI*DMOD(FPB(J)*T,1.0D0) 

XPB - XPB + SIN(WT) 

ENDDO 

DO J - 1,NSB2 

WT - TWOPI*DMOD(FSB2(J)*T,1.0D0) 

XSB2 - XSB2 + SIN(WT) 

ENDDO 

DTOT(I) - XSB1 + XPB + XSB2 
DPB(I) - XPB 
ENDDO 

! READ FILTER COEFFICIENTS 
CALL READH(NFILT,H,FILENAME,1) 

! INITIALIZE TIME DOMAIN FILTER AND PERFORM FILTERING 
DFILT(l) - FILT1(DTOT(1) ,XSAV,H,NFILT,NPT,0) 

CPUTD - SECNDS(0.0) 

DO I - 1,NDATA 

DFILT(I) - FILT1(DTOT(I),XSAV,H,NFILT,NPT,1) 

ENDDO 

CPUTD - SECNDS(CPUTD) 

! TAKE FFT OF FILTER IMPULSE RESPONSE TO GET FILTER SPECTRUM 
! AND CARRY OUT ADDITIONAL PREPARATIONS REQUIRED FOR PERFORMING 
! FREQUENCY DOMAIN FILTERING. 

DO I - 1,NFILT ! Copy the filter coefficients 

F (I) - H (I) 

ENDDO 
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e-* tt> th 


DO I * NFILT+1,NH 
F (I) - 0.0 
ENDDO 

DO I - NH+1,NHT2 
F (I) = 0.0 
ENDDO 

CALL FAST(F,NHT2) 

DO I * NH+1, NHT2 
R (I) - 0.0 
ENDDO 

! PERFORM FREQUENCY DOMAIN FILTERING BY LOOPING OVER 
! SUCCESSIVE BLOCKS OF DATA 

IF (JMOD(NDATA,NH).EQ.0) THEN 
NBLOCKS - NDATA/NH 
IPART = 0 
ELSE 

NBLOCKS - NDATA/NH + 1 
IPART - 1 
ENDIF 

CPUFD = SECNDS(0.0) 


! Pad out filter coefficients with zero 

! Force second half of filter to zero 

! Get filter spectrum 

! Zero second half of R() prior to use 


DO IB - 1,NBLOCKS 
IPOS » (IB-1)*NH 
! Transfer data for each block - 
DO I * 1,NH 

R(I) = DTOT(IPOS+I) 

ENDDO 

! Filter block of data 
CALL RFILT(R,F,WORKl,NH) 

! Save filtered data 
DO I = 1,NH 

DFIL(IPOS+I) - R(I) 

ENDDO 

ENDDO 

CPUFD - SECNDS(CPUFD) 

! WRITE OUT RESULTS TO A FILE 


WRITE(* 1 (//IX A /) *) 

& ’WRITING FILTERED SEQUENCES TO FILE TESTFILT.OUT’ 


OPEN(UNIT-1,FILE''TESTFILT.OUT’,STATUS" 
& CARRIAGECONTROL-'LIST') 


'UNKNOWN' 


) 'FILTER COEFFICIENT FILE FILENAME 

) 'NO. OF "FILTER COEFFICIENTS :',NFILT 
) 'FREQUENCIES IN 1st STOPBAND:’, 

(FSB1(I),I=1,NSB1) 

) 'FREQUENCIES IN PASSBAND :', 

(FPB(I),I-1,NPB> 

) 'FREQUENCIES IN 2nd STOPBAND:', 

(FSB2(I),I«1,NSB2) 

(2 (/, IX, A, F10.3, A) ) ') 

'CPU TIME FOR FREQ DOMAIN FILTERING - ',CPUFD,' secs', 

'CPU TIME FOR TIME DOMAIN FILTERING - ’,CPUTD,’ secs’ 
WRITEU, ' (/, IX,A,/, IX,A) ') 

'INDEX TIME X ORIG X PB X FIL TIME X FIL FREQ', 


WRITEU, 
WRITEU, 
WRITE (1, 

; 

WRITEU, 

WRITEU, 

i 

WRITE(1, 


(/, IX, A, IX, A) 
(/, IX, A, 15, /) 
(IX,A,6F7.1) ' 

(IX, A, 6F7.1) ' 

(IX,A,6F7.1) ' 


DO I - 1,NDATA 

WRITEU, ' < IX, 15, 3F11.4,2F11.5) ’) 

& I,I*DT,dTOT(I),DPB(I),DFIL(I),DFILT(I) 

ENDDO 

CLOSE(UNIT-1) 

STOP 

END 
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APPENDIX E 


A Fortran Program For Studying The Effects of FIR 
Bandpass Filtering on Pulse-Type Transients 


The following Fortran program can be used to study the effects of FIR bandpass filtering on a 
simulated pulse-like transient that is typical of vibration responses on the F/A-18. 

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

PROGRAM CHECKQ 
IMPLICIT NONE 

INTEGER NFFT,NFFT2,NPTS,NFAST 
REAL SRATE,PI,TWOPI,DT 


PARAMETER ( NFFT = 2048 ) 
PARAMETER ( NFFT2 = 2 *NFFT ) 
PARAMETER ( NFAST * NFFT*2+2 ) 
PARAMETER ( NPTS = 3000 ) 
PARAMETER ( SRATE - 606.06 ) 
PARAMETER ( DT = 1.0/SRATE } 
PARAMETER ( PI = 3.141592654 ) 
PARAMETER ( TWOPI = 2.0*PI ) 


REAL XSAV(NFFT2),H(NFFT) 

REAL XDATA(NPTS),XFILT(NPTS),WINDOW(NPTS) 

REAL FFTFILT(NFAST),FFTDATA(NFAST) 

REAL FSB1(10),FPB(10),FSB2(10) 

REAL XSB1,XSB2,XPB,OMEGA,GAMMA 

REAL MagFac,PSDfac,NoiseBW,MFFTDATA,MFFTFILT,FREQ,T 

INTEGER NFILT,INIT,IFILTER,I,J,N,NSB1,NSB2,NPB,NWIN 
INTEGER WindowCode,IPOS,IOUT 

CHARACTER ASCII1*4,ASCII2*15,ASCII3M 
CHARACTER FILENAME*20,CHOICE*20,WindowDes C*80 

REAL FILT1,COS,SIN,FLOAT,CABS,EXP 

INTEGER LENSTR 
COMPLEX CMPLX 

IFILTER = 0 

DO WHILE (IFILTER.EQ.0) 


WRITE (*,*) 

WRITE(*,*) 'PROGRAM TO TEST TRANSIENT RESPONSE OF BANDPASS '// 

5 'FILTER DESIGNS’ 

WRITE (*,*) '-'// 

6 i-. 

WRITE)*,*) ' 1 = 10 HZ TO 20 HZ BANDPASS FILTER’ 

WRITE(*,*) ' 2 » 32 HZ TO 52 HZ BANDPASS FILTER' 

WRITE!*,*) 


WRITE(*,'(1X,A,$)') 'INPUT YOUR SELECTION: ' 

READ(*,*) IFILTER 

IF (IFILTER.NE.1 .AND. IFILTER.NE.2) IFILTER - 0 
ENDDO 

WindowCode - 0 

DO WHILE (WindowCode.EQ.0) 

WRITE(*,*) 

WRITE(*,*) 'SELECT WINDOW TO BE APPLIED TO TRANSIENT' 

WRITE!*,*) '-- 

WRITE(*,*) ' 1 - SINE SQUARED WINDOW' 

WRITE(*,*) ' 2 - EXPONENTIAL WINDOW’ 

WRITE(*,*) ' 3 = UNIFORM WINDOW’ 

WRITE(*,*) 

WRITE(*,’(IX,A,$)') 'INPUT YOUR SELECTION: ' 

READ < *,*) WindowCode 
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IF (WindowCode.LT.1 .OR. WindowCode.GT.3) WindowCode - 0 
ENDDO 

IF (IFILTER.EQ.1) THEN 
FILENAME - 'BP15.0UT' 

NSB1 - 0 

NPB - 1 

NSB2 - 0 

FPB(l) - 15.0 

NWIN - 627 

ELSE IF tIFILTER.EQ.2) THEN 
FILENAME - 'BP45.0UT' 

NSB1 - 0 

NPB - 1 

NSB2 - 0 

FPB(l) * 45.0 

NWIN - 613 

END IF 

! READ FILTER COEFFICIENTS 
CALL READH(NFILT,H,FILENAME,1) 

! OPEN FILE FOR WRITING FILTER OUTPUT 
WRITE(*, ' (/, IX,A,/) ') 

& 'WRITING RESULTS TO FILE CHECKQ.OUT' 


IOUT - 1 

OPEN(UNIT-IOUT,FILE-’CHECKQ.OUT',STATUS-'UNKNOWN', 
& CARRIAGECONTROL-'LIST') 


WRITE(IOUT,'(/,IX,2A,')' ) 
WRITE : T OUT, ' (IX,A,6F7.1) ') 

WRITE(IOUT,'(IX,A,6F7.1)') 

WRITE(IOUT,'(IX,A,6F7.1)') 

WRITE(IOUT,*) 


'FILTER FILE IS ',FILENAME 
'FREQUENCIES IN 1st STOPBAND:', 
(FSB1(I),I—1, NSB1) 

•FREQUENCIES IN PASSBAND :', 
(FPB(I),I-1,NPB) 

■FREOTTENCIES IN 2nd STOPBAND:’, 
C _,bciI),I«x,NSB2) 


! Set up data points for a cosine wave, ensuring that the unity 
! value occurs at the NWIN/2+1 point, which is the centre of the 
! transient that we are trying to simulate. 

j 

! NWIN must be an odd number. 


DO I - 1,NWIN 

T - FLOAT(I-(NWIN/2+1))*DT 

XSB1 - 0.0 

XPB -0.0 

XSB2 - 0.0 

DO J - 1,NSB1 

XSB1 - XSB1 + COS(TWOPI*FSBl(J)*T) 

ENDDO 

DO J - 1,NPB 

XPB - XPB + COS(TWOPI*FPB(J)*T) 

ENDDO 

DO J - 1,NSB2 

XSB2 - XSB2 + COS(TWOPI*FSB2(J)*T) 

ENDDO 

XDATA(I ) - XSB1 + XPB + XSB2 
ENDDO 

DO I - NWIN+1,NPTS 
XDATA(I) - 0.0 
ENDDO 

IF (WindowCode.EQ.l) THEN 

! Set up a sine squared window which is similar to a 
! Hanning window, except that both the end points are 
! zero. It is created by squaring the values for the 
! first half cycle of a sine wave. 

DO I - 1,NWIN 

WINDOW(I) - ( SIN(FLOAT(I—1)*Pi/FLOAT(NWIN-1)) )**2 
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ENDDO 

WindowDesc 'COSINE SEQUENCE WITH SINE SQUARED WINDOW' 

ELSE IF (WindowCode.EQ.2) THEN 

! Set up an exponential window, increasing up to the centre of 
! the window at NWIN/2+1, and decaying from then onwards. The 
! damping is specified by GAMMA as % of critical. 

GAMMA - 0.03 
OMEGA - TWOPI*FPB(1) 

DO I - 1,NWIN 

T - FLOAT(I-(NWIN/2+1))*DT 
IF (I.LE.NWIN/2+1) THEN 

WINDOW(I) = EXP(+GAMMA*OMEGA*T) 

ELSE 

WINDOW(I) - EXP(-GAMMA*OMEGA*T) 

END IF 
ENDDO 

WindowDesc » 'COSINE SEQUENCE WITH EXPONENTIAL WINDOW' 

ELSE 

! Set up uniform window 
DO I * 1,NWIN 
WINDOW(I) - 1.0 
ENDDO 

WindowDesc * 'COSINE SEQUENCE WITH UNIFORM WINDOW’ 

END IF 

! Apply the window to the NWIN points in the data sequence 
DO I * 1,NWIN 

XDATA(I) - XDATA(I)*WINDOW(I) 

ENDDO 

! Initialize filter 

XFILT(l) = FILT1(XDATA(1),XSAV,H,NFILT,IPOS,0) 

! Filter the sequence 
DO I = 1,NPTS 

XFILT(I) - FILT1(XDATA(I),XSAV,H,NFILT,IPOS,1) 

ENDDO 

! Get FFT of raw and filtered data for comparison 

DO I = 1,NFFT 

FFTDATA(I) - XDATA(I) 

FFTFILT(I) - XFILT(I) 

ENDDO 


CALL FAST(FFTDATA,NFFT) 
CALL FAST(FFTFILT,NFFT) 


WRITE(IOUT,*) 

WRITE(IOUT,*) 

WRITE(IOUT,*) 

WRITE(IOUT,*) 

WRITE (IOUT, *) 

WRITE(IOUT,*) 

WRITE (IOUT, ' (IX, 14, A) ') 
WRITE(IOUT,*) 

WRITE(IOUT,*) 

WRITE(IOUT,*) 

WRITE(IOUT,*) 


'FAST FOURIER TRANSFORM OF RAW AND FILTERED DATA' 
' FOR CHECKING Q OF FIR FILTER' 


NFFT,' POINT FAST FOURIER TRANSFORM' 
WindowDesc(1:LENSTR(WindowDesc)) 


WRITE(IOUT,*) 'POINT FREQUENCY 


UNFILTERED 
UNFILTERED 
FFT MAG 
PSD M A 2/Hz 


FILTERED 
FILTERED 
FFT MAG 
PSD M / '2/Hz 


WRITE(IOUT,*) 


// 

// 

7/ 


DO I * 1,NFFT/2+1 
J - 2*1 

FREQ « FLOAT(1-1)/FLOAT(NFFT)*SRATE 

MFFTDATA - ( CABS( CMPLX(FFTDATA(J-l),FFTDATA(J)) ) ) 
MFFTFILT » ( CABS( CMPLX(FFTFILT(J-l),FFTFILT(J)) ) ) 
! Display FFT results as well as PSD results 
WRITE(IOUT,100) I,FREQ,MFFTDATA,MFFTFILT, 
i MFFTDATA**2/SRATE,MFFTFILT**2/SRATE 
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ENDDO 


100 FORMAT(IX,I5,F13.4,4(IX,1PE11.4)) 

WRITE(IOUT,*) 

WRITE (IOUT, *) ' ———————————————— 

WRITE(IOUT,*) ' RAW AND FILTERED TIME HISTORY DATA USED' 
WRITE(IOUT,*) ' FOR CHECKING Q OF FIR FILTER' 

WRITE (IOUT, *) ' —————————————— 

WRITE(IOUT,*) 

WRITE(IOUT,'(IX,14,A)') NFILT,' POINT BAND PASS FILTER’ 
WRITE(IOUT,*) 


WRITE(IOUT,*) 

& 

I 

' FILTERED' 



'// 

WRITE(IOUT,*) 

& 

1 

* SEQUENCE* 

UNFILTERED 

FILTERED 

'll 

WRITE(IOUT,*) 

& 

WRITE(IOUT,*) 

'POINT TIME 

' - DELAY' 

SEQUENCE 

SEQUENCE 

'll 

'll 


& 


DO I = 1,NWIN+NFILT 
T - FLOAT(1-1)/SRATE 

WRITE(IOUT,101) I,T,XDATA(I),XFILT(I),XFILT(1+(NFILT-1)/2) 
ENDDO 

101 FORMAT(IX,15,4F12.5) 

CLOSE(UNIT-IOUT) 

STOP 

END 
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